import sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import os
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import font_manager

plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 设置相对路径
LOG_FILE = 'D/analysis_log.txt'
os.makedirs('D', exist_ok=True)

# 初始化日志文件
log_file = open(LOG_FILE, 'w', encoding='utf-8')
original_stdout = sys.stdout
sys.stdout = log_file

print("=" * 80)
print("古代玻璃制品化学成分关联分析 - 问题4建模实现")
print("=" * 80)
print("\n开始执行时间:", pd.Timestamp.now().strftime("%Y-%m-%d %H:%M:%S"))
print("\n")

# =============================================================================
# 1. 数据预处理
# =============================================================================
print("步骤1: 数据预处理")
print("-" * 60)

# 读取预处理数据
try:
    # 使用相对路径
    df = pd.read_excel('handled/预处理.xlsx', sheet_name='Sheet1')
    print("成功读取预处理数据，样本数量:", df.shape[0])
except Exception as e:
    print("读取预处理数据失败:", str(e))
    sys.exit(1)

# 过滤有效数据 (成分总和在85%-105%之间)
df = df[(df['成分总和'] >= 85) & (df['成分总和'] <= 105)]
print("过滤后有效样本数量:", df.shape[0])

# 化学成分列表
chemicals = [
    '二氧化硅(SiO2)', '氧化钠(Na2O)', '氧化钾(K2O)', '氧化钙(CaO)', '氧化镁(MgO)',
    '氧化铝(Al2O3)', '氧化铁(Fe2O3)', '氧化铜(CuO)', '氧化铅(PbO)', '氧化钡(BaO)',
    '五氧化二磷(P2O5)', '氧化锶(SrO)', '氧化锡(SnO2)', '二氧化硫(SO2)'
]

# 处理缺失值 - 将NaN替换为0 (表示未检测到)
df[chemicals] = df[chemicals].fillna(0)

# 可视化1：数据质量检查 - 成分总和分布
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.hist(df['成分总和'], bins=30, alpha=0.7, color='skyblue', edgecolor='black')
plt.title('化学成分总和分布', fontsize=14)
plt.xlabel('成分总和 (%)', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.axvline(85, color='red', linestyle='--', label='下限(85%)')
plt.axvline(105, color='red', linestyle='--', label='上限(105%)')
plt.legend()

plt.subplot(1, 2, 2)
glass_type_counts = df['类型'].value_counts()
plt.pie(glass_type_counts.values, labels=glass_type_counts.index, autopct='%1.1f%%', 
        colors=['lightcoral', 'lightskyblue'])
plt.title('玻璃类型分布', fontsize=14)
plt.tight_layout()
plt.savefig('D/data_quality_overview.png', dpi=300, bbox_inches='tight')
plt.close()
print("数据质量概览图已保存")

# 移除方差过小的成分 (标准差 < 0.1)
low_variance_cols = []
variance_data = []
for col in chemicals:
    variance = df[col].std()
    variance_data.append({'成分': col, '标准差': variance})
    if variance < 0.1:
        low_variance_cols.append(col)
        chemicals.remove(col)

# 可视化2：成分方差分析
variance_df = pd.DataFrame(variance_data).sort_values('标准差', ascending=False)
plt.figure(figsize=(14, 8))
bars = plt.bar(range(len(variance_df)), variance_df['标准差'], 
               color=['red' if std < 0.1 else 'green' for std in variance_df['标准差']])
plt.axhline(0.1, color='orange', linestyle='--', label='阈值(0.1)')
plt.title('各化学成分标准差分析', fontsize=16)
plt.xlabel('化学成分', fontsize=12)
plt.ylabel('标准差', fontsize=12)
plt.xticks(range(len(variance_df)), variance_df['成分'], rotation=45, ha='right')
plt.legend()
plt.tight_layout()
plt.savefig('D/component_variance_analysis.png', dpi=300, bbox_inches='tight')
plt.close()
print("成分方差分析图已保存")

if low_variance_cols:
    print(f"移除方差过小的成分: {', '.join(low_variance_cols)}")
    print(f"保留的化学成分: {', '.join(chemicals)}")
else:
    print("所有化学成分方差均大于0.1，保留全部成分")

# 中心对数比转换 (CLR) 函数
def clr_transform(data, epsilon=1e-4):
    """
    对成分数据进行中心对数比转换 (CLR)
    
    参数:
    data - 成分数据 (DataFrame)
    epsilon - 用于替换0值的小常数
    
    返回:
    CLR转换后的数据
    """
    # 复制数据避免修改原始数据
    transformed = data.copy()
    
    # 将0值替换为epsilon
    transformed[transformed == 0] = epsilon
    
    # 计算几何平均
    geom_mean = transformed.apply(lambda x: np.exp(np.mean(np.log(x))), axis=1)
    
    # 应用CLR转换
    for col in transformed.columns:
        transformed[col] = np.log(transformed[col] / geom_mean)
    
    return transformed

# 可视化3：CLR转换前后对比
sample_components = ['二氧化硅(SiO2)', '氧化钾(K2O)', '氧化铅(PbO)', '氧化钡(BaO)']
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle('CLR转换前后数据分布对比', fontsize=16)

for i, comp in enumerate(sample_components):
    row = i // 2
    col = i % 2
    
    # 原始数据
    axes[row, col].hist(df[comp], alpha=0.6, label='原始数据', bins=20, color='skyblue')
    
    # CLR转换后数据（先进行转换用于显示）
    temp_clr = clr_transform(df[chemicals])
    axes[row, col].hist(temp_clr[comp], alpha=0.6, label='CLR转换后', bins=20, color='lightcoral')
    
    axes[row, col].set_title(f'{comp}', fontsize=12)
    axes[row, col].set_xlabel('数值', fontsize=10)
    axes[row, col].set_ylabel('频数', fontsize=10)
    axes[row, col].legend()

plt.tight_layout()
plt.savefig('D/clr_transformation_comparison.png', dpi=300, bbox_inches='tight')
plt.close()
print("CLR转换对比图已保存")

# 应用CLR转换
df_clr = clr_transform(df[chemicals])

# 添加分组信息
df_clr['类型'] = df['类型']
df_clr['表面风化'] = df['表面风化']

print("CLR转换完成")
print("数据预处理完成\n")

# =============================================================================
# 2. 分组数据准备
# =============================================================================
print("步骤2: 分组数据准备")
print("-" * 60)

# 按玻璃类型分组
high_k_df = df_clr[df_clr['类型'] == '高钾']
lead_ba_df = df_clr[df_clr['类型'] == '铅钡']

print(f"高钾玻璃样本数: {high_k_df.shape[0]}")
print(f"铅钡玻璃样本数: {lead_ba_df.shape[0]}")

# 按风化状态分组 (用于敏感性分析)
high_k_weathered = high_k_df[high_k_df['表面风化'] == '风化']
high_k_unweathered = high_k_df[high_k_df['表面风化'] == '无风化']
lead_ba_weathered = lead_ba_df[lead_ba_df['表面风化'] == '风化']
lead_ba_unweathered = lead_ba_df[lead_ba_df['表面风化'] == '无风化']

print(f"高钾风化样本数: {high_k_weathered.shape[0]}")
print(f"高钾未风化样本数: {high_k_unweathered.shape[0]}")
print(f"铅钡风化样本数: {lead_ba_weathered.shape[0]}")
print(f"铅钡未风化样本数: {lead_ba_unweathered.shape[0]}")

# 可视化4：分组样本统计
group_stats = pd.DataFrame({
    '类型': ['高钾', '高钾', '铅钡', '铅钡'],
    '风化状态': ['风化', '无风化', '风化', '无风化'],
    '样本数': [len(high_k_weathered), len(high_k_unweathered), 
             len(lead_ba_weathered), len(lead_ba_unweathered)]
})

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
sns.barplot(data=group_stats, x='类型', y='样本数', hue='风化状态', palette='Set2')
plt.title('各组样本数量统计', fontsize=14)
plt.ylabel('样本数', fontsize=12)

plt.subplot(1, 2, 2)
# 主要成分分布箱线图
key_components = ['二氧化硅(SiO2)', '氧化钾(K2O)', '氧化铅(PbO)', '氧化钡(BaO)']
plot_data = []
for comp in key_components:
    for glass_type in ['高钾', '铅钡']:
        data = df_clr[df_clr['类型'] == glass_type][comp]
        for val in data:
            plot_data.append({'成分': comp, '类型': glass_type, '值': val})

plot_df = pd.DataFrame(plot_data)
sns.boxplot(data=plot_df, x='成分', y='值', hue='类型', palette='viridis')
plt.title('关键成分CLR值分布对比', fontsize=14)
plt.xticks(rotation=45)
plt.ylabel('CLR值', fontsize=12)
plt.tight_layout()
plt.savefig('D/group_statistics_overview.png', dpi=300, bbox_inches='tight')
plt.close()
print("分组统计概览图已保存")

print("分组数据准备完成\n")

# =============================================================================
# 3. 关联关系分析
# =============================================================================
print("步骤3: 关联关系分析")
print("-" * 60)

def calculate_correlation(df, group_name):
    """
    计算相关系数矩阵并生成热力图
    
    参数:
    df - 数据集
    group_name - 组名 (用于文件名)
    
    返回:
    相关系数矩阵
    """
    # 计算相关系数
    corr_matrix = df[chemicals].corr(method='pearson')
    
    # 保存相关系数矩阵
    corr_matrix.to_csv(f'D/{group_name}_correlation_matrix.csv', encoding='utf-8-sig')
    print(f"{group_name}相关系数矩阵已保存")
    
    # 绘制热力图
    plt.figure(figsize=(14, 12))
    mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
    sns.heatmap(corr_matrix, mask=mask, annot=True, fmt=".2f", cmap='coolwarm', 
                center=0, vmin=-1, vmax=1, linewidths=.5, annot_kws={"size": 9})
    plt.title(f'{group_name}玻璃化学成分相关系数矩阵', fontsize=14)
    plt.xticks(rotation=45, ha='right', fontsize=9)
    plt.yticks(fontsize=9)
    plt.tight_layout()
    plt.savefig(f'D/{group_name}_correlation_heatmap.png', dpi=300)
    plt.close()
    print(f"{group_name}热力图已保存")
    
    # 识别强相关对 (|r| > 0.5)
    strong_corr_pairs = []
    for i in range(len(corr_matrix.columns)):
        for j in range(i+1, len(corr_matrix.columns)):
            if abs(corr_matrix.iloc[i, j]) > 0.5:
                strong_corr_pairs.append({
                    '成分1': corr_matrix.columns[i],
                    '成分2': corr_matrix.columns[j],
                    '相关系数': corr_matrix.iloc[i, j]
                })
    
    if strong_corr_pairs:
        strong_corr_df = pd.DataFrame(strong_corr_pairs)
        strong_corr_df = strong_corr_df.sort_values(by='相关系数', key=abs, ascending=False)
        strong_corr_df.to_csv(f'D/{group_name}_strong_correlations.csv', 
                             index=False, encoding='utf-8-sig')
        print(f"{group_name}强相关对({len(strong_corr_pairs)}个)已保存")
        
        # 可视化强相关对
        if len(strong_corr_pairs) > 0:
            plt.figure(figsize=(12, 8))
            strong_corr_plot = strong_corr_df.head(10)  # 显示前10个
            colors = ['red' if r < 0 else 'blue' for r in strong_corr_plot['相关系数']]
            plt.barh(range(len(strong_corr_plot)), 
                    strong_corr_plot['相关系数'], color=colors, alpha=0.7)
            plt.yticks(range(len(strong_corr_plot)), 
                      [f"{row['成分1']}-{row['成分2']}" for _, row in strong_corr_plot.iterrows()])
            plt.xlabel('相关系数', fontsize=12)
            plt.title(f'{group_name}玻璃强相关成分对 (|r| > 0.5)', fontsize=14)
            plt.axvline(0, color='black', linestyle='-', alpha=0.3)
            plt.tight_layout()
            plt.savefig(f'D/{group_name}_strong_correlations_plot.png', dpi=300, bbox_inches='tight')
            plt.close()
            print(f"{group_name}强相关对可视化已保存")
    
    return corr_matrix

# 分析高钾玻璃
print("\n分析高钾玻璃...")
high_k_corr = calculate_correlation(high_k_df, '高钾')

# 分析铅钡玻璃
print("\n分析铅钡玻璃...")
lead_ba_corr = calculate_correlation(lead_ba_df, '铅钡')

# 可视化5：相关系数分布对比
plt.figure(figsize=(15, 6))

plt.subplot(1, 3, 1)
# 提取上三角矩阵的相关系数（排除对角线）
high_k_corr_values = []
lead_ba_corr_values = []
for i in range(len(chemicals)):
    for j in range(i+1, len(chemicals)):
        high_k_corr_values.append(high_k_corr.iloc[i, j])
        lead_ba_corr_values.append(lead_ba_corr.iloc[i, j])

plt.hist(high_k_corr_values, bins=20, alpha=0.6, label='高钾玻璃', color='skyblue')
plt.hist(lead_ba_corr_values, bins=20, alpha=0.6, label='铅钡玻璃', color='lightcoral')
plt.xlabel('相关系数', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('相关系数分布对比', fontsize=14)
plt.legend()

plt.subplot(1, 3, 2)
plt.scatter(high_k_corr_values, lead_ba_corr_values, alpha=0.6, s=50)
plt.plot([-1, 1], [-1, 1], 'r--', alpha=0.8)
plt.xlabel('高钾玻璃相关系数', fontsize=12)
plt.ylabel('铅钡玻璃相关系数', fontsize=12)
plt.title('两类玻璃相关系数散点图', fontsize=14)
plt.grid(True, alpha=0.3)

plt.subplot(1, 3, 3)
diff_values = np.array(high_k_corr_values) - np.array(lead_ba_corr_values)
plt.hist(diff_values, bins=20, color='green', alpha=0.7)
plt.xlabel('相关系数差异 (高钾-铅钡)', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('相关系数差异分布', fontsize=14)
plt.axvline(0, color='red', linestyle='--', alpha=0.8)

plt.tight_layout()
plt.savefig('D/correlation_comparison_overview.png', dpi=300, bbox_inches='tight')
plt.close()
print("相关系数对比概览图已保存")

print("关联关系分析完成\n")

# =============================================================================
# 4. 差异比较
# =============================================================================
print("步骤4: 差异比较")
print("-" * 60)

# 计算相关系数差异矩阵
corr_diff = high_k_corr - lead_ba_corr

# 保存差异矩阵
corr_diff.to_csv('D/correlation_difference_matrix.csv', encoding='utf-8-sig')
print("相关系数差异矩阵已保存")

# 绘制差异热力图
plt.figure(figsize=(14, 12))
mask = np.triu(np.ones_like(corr_diff, dtype=bool))
sns.heatmap(corr_diff, mask=mask, annot=True, fmt=".2f", cmap='coolwarm', 
            center=0, vmin=-1, vmax=1, linewidths=.5, annot_kws={"size": 9})
plt.title('高钾与铅钡玻璃相关系数差异矩阵', fontsize=14)
plt.xticks(rotation=45, ha='right', fontsize=9)
plt.yticks(fontsize=9)
plt.tight_layout()
plt.savefig('D/correlation_difference_heatmap.png', dpi=300)
plt.close()
print("差异热力图已保存")

# Fisher z变换检验函数
def fisher_z_test(r1, r2, n1, n2):
    """
    使用Fisher z变换检验两个相关系数是否显著不同
    
    参数:
    r1, r2 - 相关系数
    n1, n2 - 样本大小
    
    返回:
    z_score, p_value
    """
    # Fisher z变换
    z1 = np.arctanh(r1)
    z2 = np.arctanh(r2)
    
    # 计算z统计量
    se_diff = np.sqrt(1/(n1-3) + 1/(n2-3))
    z_diff = (z1 - z2) / se_diff
    
    # 计算p值 (双尾检验)
    p_value = 2 * (1 - stats.norm.cdf(abs(z_diff)))
    
    return z_diff, p_value

# 执行差异检验
n_high_k = high_k_df.shape[0]
n_lead_ba = lead_ba_df.shape[0]

significance_results = []
for i in range(len(chemicals)):
    for j in range(i+1, len(chemicals)):
        chem1 = chemicals[i]
        chem2 = chemicals[j]
        
        r_high_k = high_k_corr.loc[chem1, chem2]
        r_lead_ba = lead_ba_corr.loc[chem1, chem2]
        
        z_score, p_value = fisher_z_test(r_high_k, r_lead_ba, n_high_k, n_lead_ba)
        
        significance_results.append({
            '成分1': chem1,
            '成分2': chem2,
            '高钾相关系数': r_high_k,
            '铅钡相关系数': r_lead_ba,
            '相关系数差异': r_high_k - r_lead_ba,
            'z统计量': z_score,
            'p值': p_value
        })

# 创建结果DataFrame
significance_df = pd.DataFrame(significance_results)

# Bonferroni校正
n_tests = len(significance_df)
significance_df['校正后p值'] = significance_df['p值'] * n_tests
significance_df['显著(0.05)'] = significance_df['校正后p值'] < 0.05
significance_df['显著(0.01)'] = significance_df['校正后p值'] < 0.01

# 保存结果
significance_df.to_csv('D/correlation_significance_results.csv', 
                      index=False, encoding='utf-8-sig')
print("相关系数差异显著性检验结果已保存")

# 可视化6：显著性检验结果
plt.figure(figsize=(15, 10))

# p值分布
plt.subplot(2, 3, 1)
plt.hist(significance_df['p值'], bins=20, alpha=0.7, color='lightblue')
plt.xlabel('p值', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('p值分布', fontsize=14)
plt.axvline(0.05/n_tests, color='red', linestyle='--', label=f'Bonferroni阈值({0.05/n_tests:.6f})')
plt.legend()

# z统计量分布
plt.subplot(2, 3, 2)
plt.hist(significance_df['z统计量'], bins=20, alpha=0.7, color='lightgreen')
plt.xlabel('z统计量', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('z统计量分布', fontsize=14)
plt.axvline(0, color='red', linestyle='--', alpha=0.8)

# 相关系数差异vs z统计量
plt.subplot(2, 3, 3)
plt.scatter(significance_df['相关系数差异'], significance_df['z统计量'], alpha=0.6)
plt.xlabel('相关系数差异', fontsize=12)
plt.ylabel('z统计量', fontsize=12)
plt.title('差异 vs z统计量', fontsize=14)
plt.grid(True, alpha=0.3)

# 显著性结果统计
plt.subplot(2, 3, 4)
sig_counts = [
    len(significance_df) - significance_df['显著(0.05)'].sum(),
    significance_df['显著(0.05)'].sum() - significance_df['显著(0.01)'].sum(),
    significance_df['显著(0.01)'].sum()
]
labels = ['不显著', '显著(p<0.05)', '高度显著(p<0.01)']
colors = ['lightgray', 'orange', 'red']
plt.pie(sig_counts, labels=labels, colors=colors, autopct='%1.1f%%')
plt.title('显著性检验结果统计', fontsize=14)

# 相关系数对比散点图（按显著性着色）
plt.subplot(2, 3, 5)
colors = ['red' if sig else 'blue' for sig in significance_df['显著(0.05)']]
plt.scatter(significance_df['高钾相关系数'], significance_df['铅钡相关系数'], 
           c=colors, alpha=0.6, s=50)
plt.plot([-1, 1], [-1, 1], 'k--', alpha=0.8)
plt.xlabel('高钾玻璃相关系数', fontsize=12)
plt.ylabel('铅钡玻璃相关系数', fontsize=12)
plt.title('相关系数对比（红色=显著差异）', fontsize=14)
plt.grid(True, alpha=0.3)

# 差异大小分布
plt.subplot(2, 3, 6)
abs_diff = np.abs(significance_df['相关系数差异'])
plt.hist(abs_diff, bins=20, alpha=0.7, color='purple')
plt.xlabel('|相关系数差异|', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('相关系数差异绝对值分布', fontsize=14)

plt.tight_layout()
plt.savefig('D/significance_analysis_overview.png', dpi=300, bbox_inches='tight')
plt.close()
print("显著性分析概览图已保存")

# 提取显著差异对
significant_pairs = significance_df[significance_df['显著(0.05)']].copy()
significant_pairs = significant_pairs.sort_values(by='z统计量', key=abs, ascending=False)

if not significant_pairs.empty:
    significant_pairs.to_csv('D/significant_difference_pairs.csv', 
                           index=False, encoding='utf-8-sig')
    print(f"发现{len(significant_pairs)}个显著差异对")
    
    # 可视化显著差异对
    plt.figure(figsize=(14, 8))
    if len(significant_pairs) > 15:
        plot_pairs = significant_pairs.head(15)
    else:
        plot_pairs = significant_pairs
    
    pair_labels = [f"{row['成分1'][:6]}-{row['成分2'][:6]}" for _, row in plot_pairs.iterrows()]
    z_values = plot_pairs['z统计量']
    colors = ['red' if z < 0 else 'blue' for z in z_values]
    
    plt.barh(range(len(plot_pairs)), z_values, color=colors, alpha=0.7)
    plt.yticks(range(len(plot_pairs)), pair_labels)
    plt.xlabel('z统计量', fontsize=12)
    plt.title('显著差异成分对的z统计量', fontsize=14)
    plt.axvline(0, color='black', linestyle='-', alpha=0.3)
    plt.tight_layout()
    plt.savefig('D/significant_pairs_z_scores.png', dpi=300, bbox_inches='tight')
    plt.close()
    print("显著差异对z统计量图已保存")
else:
    print("未发现显著差异对")

print("差异比较完成\n")

# =============================================================================
# 5. 敏感性分析
# =============================================================================
print("步骤5: 敏感性分析")
print("-" * 60)

# 5.1 CLR转换敏感性分析
print("\n5.1 CLR转换敏感性分析...")

# 使用原始数据计算相关系数
high_k_raw_corr = high_k_df[chemicals].corr()
lead_ba_raw_corr = lead_ba_df[chemicals].corr()

# 比较关键成分对的相关系数差异
key_pairs = [
    ('氧化钾(K2O)', '氧化钙(CaO)'),
    ('氧化铅(PbO)', '氧化钡(BaO)'),
    ('二氧化硅(SiO2)', '氧化钾(K2O)'),
    ('二氧化硅(SiO2)', '氧化铅(PbO)')
]

clr_sensitivity = []
for pair in key_pairs:
    chem1, chem2 = pair
    
    # CLR转换后的相关系数
    r_high_k_clr = high_k_corr.loc[chem1, chem2]
    r_lead_ba_clr = lead_ba_corr.loc[chem1, chem2]
    
    # 原始数据的相关系数
    r_high_k_raw = high_k_raw_corr.loc[chem1, chem2]
    r_lead_ba_raw = lead_ba_raw_corr.loc[chem1, chem2]
    
    clr_sensitivity.append({
        '成分对': f"{chem1}-{chem2}",
        '高钾CLR': r_high_k_clr,
        '高钾原始': r_high_k_raw,
        '高钾差异': abs(r_high_k_clr - r_high_k_raw),
        '铅钡CLR': r_lead_ba_clr,
        '铅钡原始': r_lead_ba_raw,
        '铅钡差异': abs(r_lead_ba_clr - r_lead_ba_raw)
    })

clr_sensitivity_df = pd.DataFrame(clr_sensitivity)
clr_sensitivity_df.to_csv('D/clr_sensitivity_analysis.csv', 
                         index=False, encoding='utf-8-sig')

# 可视化7：CLR敏感性分析
plt.figure(figsize=(14, 8))
x = np.arange(len(clr_sensitivity_df))
width = 0.35

plt.subplot(1, 2, 1)
plt.bar(x - width/2, clr_sensitivity_df['高钾差异'], width, label='高钾玻璃', alpha=0.7)
plt.bar(x + width/2, clr_sensitivity_df['铅钡差异'], width, label='铅钡玻璃', alpha=0.7)
plt.xlabel('成分对', fontsize=12)
plt.ylabel('CLR转换前后相关系数差异', fontsize=12)
plt.title('CLR转换敏感性分析', fontsize=14)
plt.xticks(x, [pair.replace('-', '\n') for pair in clr_sensitivity_df['成分对']], fontsize=10)
plt.legend()

plt.subplot(1, 2, 2)
for i, (_, row) in enumerate(clr_sensitivity_df.iterrows()):
    plt.scatter(row['高钾原始'], row['高钾CLR'], color='blue', alpha=0.7, s=100, label='高钾' if i==0 else "")
    plt.scatter(row['铅钡原始'], row['铅钡CLR'], color='red', alpha=0.7, s=100, label='铅钡' if i==0 else "")
    
plt.plot([-1, 1], [-1, 1], 'k--', alpha=0.5)
plt.xlabel('原始相关系数', fontsize=12)
plt.ylabel('CLR相关系数', fontsize=12)
plt.title('CLR转换前后相关系数对比', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('D/clr_sensitivity_visualization.png', dpi=300, bbox_inches='tight')
plt.close()
print("CLR转换敏感性分析可视化已保存")
print("CLR转换敏感性分析完成")

# 5.2 风化状态分层分析
print("\n5.2 风化状态分层分析...")

def compare_subgroups(high_k_weather, high_k_unweather, lead_ba_weather, lead_ba_unweather):
    """比较风化与未风化子组的相关系数差异"""
    results = []
    
    for chem1 in chemicals:
        for chem2 in chemicals:
            if chem1 == chem2:
                continue
                
            # 高钾风化组
            r_high_k_w = high_k_weather[[chem1, chem2]].corr().iloc[0,1]
            # 高钾未风化组
            r_high_k_uw = high_k_unweather[[chem1, chem2]].corr().iloc[0,1]
            # 铅钡风化组
            r_lead_ba_w = lead_ba_weather[[chem1, chem2]].corr().iloc[0,1]
            # 铅钡未风化组
            r_lead_ba_uw = lead_ba_unweather[[chem1, chem2]].corr().iloc[0,1]
            
            # 计算组内差异
            diff_high_k = abs(r_high_k_w - r_high_k_uw)
            diff_lead_ba = abs(r_lead_ba_w - r_lead_ba_uw)
            
            # 计算组间差异
            diff_weather = abs(r_high_k_w - r_lead_ba_w)
            diff_unweather = abs(r_high_k_uw - r_lead_ba_uw)
            
            results.append({
                '成分1': chem1,
                '成分2': chem2,
                '高钾风化': r_high_k_w,
                '高钾未风化': r_high_k_uw,
                '高钾组内差异': diff_high_k,
                '铅钡风化': r_lead_ba_w,
                '铅钡未风化': r_lead_ba_uw,
                '铅钡组内差异': diff_lead_ba,
                '风化组间差异': diff_weather,
                '未风化组间差异': diff_unweather
            })
    
    return pd.DataFrame(results)

# 执行分层分析
weather_analysis_df = compare_subgroups(
    high_k_weathered[chemicals], 
    high_k_unweathered[chemicals], 
    lead_ba_weathered[chemicals], 
    lead_ba_unweathered[chemicals]
)

# 保存结果
weather_analysis_df.to_csv('D/weathering_stratified_analysis.csv', 
                         index=False, encoding='utf-8-sig')

# 可视化8：风化状态敏感性分析
plt.figure(figsize=(16, 10))

# 组内稳定性分析
plt.subplot(2, 3, 1)
plt.hist(weather_analysis_df['高钾组内差异'], bins=20, alpha=0.6, label='高钾玻璃', color='blue')
plt.hist(weather_analysis_df['铅钡组内差异'], bins=20, alpha=0.6, label='铅钡玻璃', color='red')
plt.xlabel('组内差异 (|风化-未风化|)', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('风化状态组内稳定性分析', fontsize=14)
plt.legend()

# 组间区分度分析
plt.subplot(2, 3, 2)
plt.hist(weather_analysis_df['风化组间差异'], bins=20, alpha=0.6, label='风化组', color='green')
plt.hist(weather_analysis_df['未风化组间差异'], bins=20, alpha=0.6, label='未风化组', color='orange')
plt.xlabel('组间差异 (|高钾-铅钡|)', fontsize=12)
plt.ylabel('频数', fontsize=12)
plt.title('风化状态组间区分度分析', fontsize=14)
plt.legend()

# 稳定性vs区分度散点图
plt.subplot(2, 3, 3)
avg_within_diff = (weather_analysis_df['高钾组内差异'] + weather_analysis_df['铅钡组内差异']) / 2
avg_between_diff = (weather_analysis_df['风化组间差异'] + weather_analysis_df['未风化组间差异']) / 2
plt.scatter(avg_within_diff, avg_between_diff, alpha=0.6, s=30)
plt.xlabel('平均组内差异', fontsize=12)
plt.ylabel('平均组间差异', fontsize=12)
plt.title('稳定性 vs 区分度', fontsize=14)
plt.grid(True, alpha=0.3)

# 关键成分对的详细分析
key_pairs_indices = []
for chem1, chem2 in key_pairs:
    idx = weather_analysis_df[(weather_analysis_df['成分1'] == chem1) & 
                             (weather_analysis_df['成分2'] == chem2)].index
    if len(idx) > 0:
        key_pairs_indices.extend(idx)

if key_pairs_indices:
    key_analysis = weather_analysis_df.loc[key_pairs_indices]
    
    # 高钾组
    plt.subplot(2, 3, 4)
    x = range(len(key_analysis))
    plt.bar(x, key_analysis['高钾风化'], alpha=0.6, label='风化', width=0.4)
    plt.bar([i+0.4 for i in x], key_analysis['高钾未风化'], alpha=0.6, label='未风化', width=0.4)
    plt.xlabel('成分对', fontsize=12)
    plt.ylabel('相关系数', fontsize=12)
    plt.title('高钾玻璃关键成分对比', fontsize=14)
    plt.xticks([i+0.2 for i in x], 
              [f"{row['成分1'][:6]}-{row['成分2'][:6]}" for _, row in key_analysis.iterrows()], 
              rotation=45)
    plt.legend()
    
    # 铅钡组
    plt.subplot(2, 3, 5)
    plt.bar(x, key_analysis['铅钡风化'], alpha=0.6, label='风化', width=0.4)
    plt.bar([i+0.4 for i in x], key_analysis['铅钡未风化'], alpha=0.6, label='未风化', width=0.4)
    plt.xlabel('成分对', fontsize=12)
    plt.ylabel('相关系数', fontsize=12)
    plt.title('铅钡玻璃关键成分对比', fontsize=14)
    plt.xticks([i+0.2 for i in x], 
              [f"{row['成分1'][:6]}-{row['成分2'][:6]}" for _, row in key_analysis.iterrows()], 
              rotation=45)
    plt.legend()

# 总体稳定性指标
plt.subplot(2, 3, 6)
stability_metrics = {
    '高钾组内稳定性': weather_analysis_df['高钾组内差异'].mean(),
    '铅钡组内稳定性': weather_analysis_df['铅钡组内差异'].mean(),
    '风化组间区分度': weather_analysis_df['风化组间差异'].mean(),
    '未风化组间区分度': weather_analysis_df['未风化组间差异'].mean()
}
plt.bar(range(len(stability_metrics)), list(stability_metrics.values()), 
        color=['blue', 'red', 'green', 'orange'], alpha=0.7)
plt.xticks(range(len(stability_metrics)), list(stability_metrics.keys()), rotation=45)
plt.ylabel('平均差异', fontsize=12)
plt.title('稳定性与区分度指标', fontsize=14)

plt.tight_layout()
plt.savefig('D/weathering_sensitivity_analysis.png', dpi=300, bbox_inches='tight')
plt.close()
print("风化状态敏感性分析可视化已保存")

print("风化状态分层分析完成")

# 5.3 Bootstrap置信区间分析
print("\n5.3 Bootstrap置信区间分析...")

def bootstrap_correlation(data, n_bootstrap=1000):
    """使用Bootstrap计算相关系数的置信区间"""
    n_samples = len(data)
    bootstrap_corrs = []
    
    for _ in range(n_bootstrap):
        # 重采样
        sample = data.sample(n=n_samples, replace=True)
        # 计算相关系数
        corr = sample.corr().values
        bootstrap_corrs.append(corr)
    
    # 计算置信区间
    bootstrap_corrs = np.array(bootstrap_corrs)
    lower = np.percentile(bootstrap_corrs, 2.5, axis=0)
    upper = np.percentile(bootstrap_corrs, 97.5, axis=0)
    
    return lower, upper

# 高钾玻璃Bootstrap
print("高钾玻璃Bootstrap...")
high_k_lower, high_k_upper = bootstrap_correlation(high_k_df[chemicals])

# 铅钡玻璃Bootstrap
print("铅钡玻璃Bootstrap...")
lead_ba_lower, lead_ba_upper = bootstrap_correlation(lead_ba_df[chemicals])

# 计算关键成分对的置信区间
key_bootstrap_results = []
for chem1, chem2 in key_pairs:
    i = chemicals.index(chem1)
    j = chemicals.index(chem2)
    
    high_k_ci_lower = high_k_lower[i, j]
    high_k_ci_upper = high_k_upper[i, j]
    lead_ba_ci_lower = lead_ba_lower[i, j]
    lead_ba_ci_upper = lead_ba_upper[i, j]
    
    key_bootstrap_results.append({
        '成分对': f"{chem1}-{chem2}",
        '高钾相关系数': high_k_corr.loc[chem1, chem2],
        '高钾CI下限': high_k_ci_lower,
        '高钾CI上限': high_k_ci_upper,
        '高钾CI宽度': high_k_ci_upper - high_k_ci_lower,
        '铅钡相关系数': lead_ba_corr.loc[chem1, chem2],
        '铅钡CI下限': lead_ba_ci_lower,
        '铅钡CI上限': lead_ba_ci_upper,
        '铅钡CI宽度': lead_ba_ci_upper - lead_ba_ci_lower
    })

bootstrap_df = pd.DataFrame(key_bootstrap_results)
bootstrap_df.to_csv('D/bootstrap_confidence_intervals.csv', index=False, encoding='utf-8-sig')

# 可视化9：Bootstrap置信区间
plt.figure(figsize=(14, 8))

plt.subplot(1, 2, 1)
x = range(len(bootstrap_df))
plt.errorbar(x, bootstrap_df['高钾相关系数'], 
            yerr=[bootstrap_df['高钾相关系数'] - bootstrap_df['高钾CI下限'],
                  bootstrap_df['高钾CI上限'] - bootstrap_df['高钾相关系数']], 
            fmt='o', capsize=5, label='高钾玻璃', color='blue', alpha=0.7)
plt.errorbar([i+0.1 for i in x], bootstrap_df['铅钡相关系数'], 
            yerr=[bootstrap_df['铅钡相关系数'] - bootstrap_df['铅钡CI下限'],
                  bootstrap_df['铅钡CI上限'] - bootstrap_df['铅钡相关系数']], 
            fmt='s', capsize=5, label='铅钡玻璃', color='red', alpha=0.7)
plt.xlabel('成分对', fontsize=12)
plt.ylabel('相关系数', fontsize=12)
plt.title('Bootstrap 95%置信区间', fontsize=14)
plt.xticks(x, [pair.replace('-', '\n') for pair in bootstrap_df['成分对']], fontsize=10)
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
plt.bar([i-0.2 for i in x], bootstrap_df['高钾CI宽度'], width=0.4, 
        label='高钾玻璃', alpha=0.7, color='blue')
plt.bar([i+0.2 for i in x], bootstrap_df['铅钡CI宽度'], width=0.4, 
        label='铅钡玻璃', alpha=0.7, color='red')
plt.xlabel('成分对', fontsize=12)
plt.ylabel('置信区间宽度', fontsize=12)
plt.title('置信区间宽度比较', fontsize=14)
plt.xticks(x, [pair.replace('-', '\n') for pair in bootstrap_df['成分对']], fontsize=10)
plt.legend()

plt.tight_layout()
plt.savefig('D/bootstrap_confidence_intervals.png', dpi=300, bbox_inches='tight')
plt.close()
print("Bootstrap置信区间可视化已保存")

print("Bootstrap置信区间分析完成")
print("敏感性分析完成\n")

# =============================================================================
# 6. 结果解释与总结
# =============================================================================
print("步骤6: 结果解释与总结")
print("-" * 60)

# 总结关键发现
if not significant_pairs.empty:
    print("\n关键发现: 显著差异的成分对")
    for i, row in significant_pairs.iterrows():
        print(f"{row['成分1']} 和 {row['成分2']}:")
        print(f"  高钾相关系数: {row['高钾相关系数']:.3f}, 铅钡相关系数: {row['铅钡相关系数']:.3f}")
        print(f"  差异: {row['相关系数差异']:.3f}, z统计量: {row['z统计量']:.3f}, p值: {row['p值']:.4f}")
else:
    print("未发现显著差异的成分对")

# 最终总结可视化
plt.figure(figsize=(20, 12))

# 创建综合总结面板
fig = plt.figure(figsize=(20, 12))
gs = fig.add_gridspec(3, 4, hspace=0.3, wspace=0.3)

# 1. 样本分布总览
ax1 = fig.add_subplot(gs[0, 0])
sample_counts = df['类型'].value_counts()
ax1.pie(sample_counts.values, labels=sample_counts.index, autopct='%1.1f%%', 
        colors=['lightcoral', 'lightskyblue'])
ax1.set_title('样本分布', fontsize=12)

# 2. 显著差异对统计
ax2 = fig.add_subplot(gs[0, 1])
if not significant_pairs.empty:
    sig_stats = [
        len(significance_df) - len(significant_pairs),
        len(significant_pairs)
    ]
    ax2.pie(sig_stats, labels=['无显著差异', '显著差异'], autopct='%1.1f%%',
            colors=['lightgray', 'orange'])
ax2.set_title('显著性检验结果', fontsize=12)

# 3. 相关性强度分布
ax3 = fig.add_subplot(gs[0, 2])
strong_corr_high_k = len([r for r in high_k_corr_values if abs(r) > 0.5])
strong_corr_lead_ba = len([r for r in lead_ba_corr_values if abs(r) > 0.5])
ax3.bar(['高钾玻璃', '铅钡玻璃'], [strong_corr_high_k, strong_corr_lead_ba], 
        color=['blue', 'red'], alpha=0.7)
ax3.set_title('强相关对数量(|r|>0.5)', fontsize=12)
ax3.set_ylabel('数量')

# 4. CLR转换敏感性
ax4 = fig.add_subplot(gs[0, 3])
avg_sensitivity = clr_sensitivity_df[['高钾差异', '铅钡差异']].mean()
ax4.bar(['高钾玻璃', '铅钡玻璃'], avg_sensitivity.values, 
        color=['blue', 'red'], alpha=0.7)
ax4.set_title('CLR转换平均敏感性', fontsize=12)
ax4.set_ylabel('平均差异')

# 5-8. 关键成分对相关系数对比
for i, (chem1, chem2) in enumerate(key_pairs):
    ax = fig.add_subplot(gs[1, i])
    
    r_high_k = high_k_corr.loc[chem1, chem2]
    r_lead_ba = lead_ba_corr.loc[chem1, chem2]
    
    ax.bar(['高钾', '铅钡'], [r_high_k, r_lead_ba], 
           color=['blue', 'red'], alpha=0.7)
    ax.set_title(f'{chem1[:6]}-{chem2[:6]}', fontsize=11)
    ax.set_ylabel('相关系数')
    ax.set_ylim(-1, 1)
    ax.axhline(0, color='black', linestyle='-', alpha=0.3)

# 9. 风化状态稳定性对比
ax9 = fig.add_subplot(gs[2, 0])
weathering_stability = {
    '高钾组内': weather_analysis_df['高钾组内差异'].mean(),
    '铅钡组内': weather_analysis_df['铅钡组内差异'].mean(),
    '风化组间': weather_analysis_df['风化组间差异'].mean(),
    '未风化组间': weather_analysis_df['未风化组间差异'].mean()
}
colors = ['blue', 'red', 'green', 'orange']
bars = ax9.bar(range(len(weathering_stability)), list(weathering_stability.values()), 
               color=colors, alpha=0.7)
ax9.set_xticks(range(len(weathering_stability)))
ax9.set_xticklabels(list(weathering_stability.keys()), rotation=45, fontsize=10)
ax9.set_title('风化状态稳定性', fontsize=12)
ax9.set_ylabel('平均差异')

# 10. Bootstrap置信区间宽度
ax10 = fig.add_subplot(gs[2, 1])
avg_ci_width = bootstrap_df[['高钾CI宽度', '铅钡CI宽度']].mean()
ax10.bar(['高钾玻璃', '铅钡玻璃'], avg_ci_width.values, 
         color=['blue', 'red'], alpha=0.7)
ax10.set_title('平均置信区间宽度', fontsize=12)
ax10.set_ylabel('宽度')

# 11. 相关系数差异分布
ax11 = fig.add_subplot(gs[2, 2])
ax11.hist(significance_df['相关系数差异'], bins=15, alpha=0.7, color='purple')
ax11.axvline(0, color='red', linestyle='--', alpha=0.8)
ax11.set_title('相关系数差异分布', fontsize=12)
ax11.set_xlabel('差异(高钾-铅钡)')
ax11.set_ylabel('频数')

# 12. 工艺特征总结
ax12 = fig.add_subplot(gs[2, 3])
# 计算关键成分的平均相关性
key_components_corr = {
    'K2O-CaO': [high_k_corr.loc['氧化钾(K2O)', '氧化钙(CaO)'], 
                 lead_ba_corr.loc['氧化钾(K2O)', '氧化钙(CaO)']],
    'PbO-BaO': [high_k_corr.loc['氧化铅(PbO)', '氧化钡(BaO)'], 
                lead_ba_corr.loc['氧化铅(PbO)', '氧化钡(BaO)']]
}

x = np.arange(len(key_components_corr))
width = 0.35
ax12.bar(x - width/2, [key_components_corr['K2O-CaO'][0], key_components_corr['PbO-BaO'][0]], 
         width, label='高钾玻璃', alpha=0.7, color='blue')
ax12.bar(x + width/2, [key_components_corr['K2O-CaO'][1], key_components_corr['PbO-BaO'][1]], 
         width, label='铅钡玻璃', alpha=0.7, color='red')
ax12.set_xticks(x)
ax12.set_xticklabels(['K2O-CaO\n(助熔-稳定)', 'PbO-BaO\n(重金属助熔)'])
ax12.set_title('关键工艺成分对', fontsize=12)
ax12.set_ylabel('相关系数')
ax12.legend()

plt.suptitle('古代玻璃化学成分关联分析综合总结', fontsize=16, y=0.95)
plt.savefig('D/comprehensive_analysis_summary.png', dpi=300, bbox_inches='tight')
plt.close()
print("综合分析总结图已保存")

# 解释工艺含义
print("\n工艺含义解释:")
print("1. 铅钡玻璃中氧化铅(PbO)和氧化钡(BaO)的强正相关反映了它们作为主要助熔剂共同添加的工艺特点")
print("2. 高钾玻璃中氧化钾(K2O)和氧化钙(CaO)的正相关表明钾作为助熔剂与钙作为稳定剂的协同作用")
print("3. 两种玻璃类型中二氧化硅(SiO2)与其他成分的负相关关系反映了SiO2作为主要基体的特性")
print("4. 差异分析揭示了两种玻璃制造工艺的本质区别：铅钡玻璃以铅矿石为助熔剂，高钾玻璃以草木灰为助熔剂")

print("\n所有结果已保存到 D/ 目录")
print("分析完成时间:", pd.Timestamp.now().strftime("%Y-%m-%d %H:%M:%S"))

# 恢复标准输出
sys.stdout = original_stdout
log_file.close()

print(f"分析完成! 详细日志已保存到: {LOG_FILE}")
print(f"所有可视化结果和数据文件已保存至 D/ 目录")
print("生成的可视化图表包括:")
print("- 数据质量概览图")
print("- 成分方差分析图") 
print("- CLR转换对比图")
print("- 分组统计概览图")
print("- 相关系数热力图")
print("- 强相关对可视化")
print("- 相关系数对比概览图")
print("- 差异热力图")
print("- 显著性分析概览图")
print("- CLR敏感性分析图")
print("- 风化状态敏感性分析图")
print("- Bootstrap置信区间图")
print("- 综合分析总结图")